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ff^ ■ Abstract 

We introduce a generalized framework for sampling and reconstruction in separable Hilbert spaces. 

■^ ' Specifically, we establish that it is always possible to stably reconstruct a vector in an arbitrary Riesz 

^^ \ basis from sufficiently many of its samples in any other Riesz basis. This framework can be viewed 

as an extension of that of Eldar et al. However, whilst the latter imposes stringent assumptions on 

the reconstruction basis, and may in practice be unstable, our framework allows for recovery in any 

2 ' (Riesz) basis in a manner that is completely stable. 

Whilst the classical Shannon Sampling Theorem is a special case of our theorem, this framework 
allows us to exploit additional information about the approximated vector (or, in this case, function), 
«vj ! for example sparsity or regularity, to design a reconstruction basis that is better suited. Examples are 

^ ' presented illustrating this procedure. 
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^ '• 1 Introduction 

The Shannon Sampling Theorem, or the Nyquist-Shannon Sampling Theorem as it is also called (we will 
refer to it as the NS-Sampling Theorem throughout the paper), is a mainstay in modem signal process- 
r> ' ing and has become one of the most important theorems in mathematics of information |[26l . The list 

C^ ■ of applications of the theorem is long, and ranges from Magnetic Resonance Imaging (MRI) to sound 

engineering. We will in this paper address the question on whether or not the NS-Sampling Theorem 
can be improved. In particular, given the same set of information, could one design a reconstruction of 
a function that would be better than that provided by the NS-Sampling Theorem. The answer to such 
a question will obviously depend on the type of functions considered. However, suppose that we have 
some extra information about the functions to be reconstructed. One may, for example, have information 
about a basis that is particularly suited for such functions. Could this information be used to improve the 
reconstruction given by the NS-Sampling Theorem, even if it is based on the same sampling procedure? 
Although such a question has been posed before, and numerous extensions of the NS-Sampling Theorem 
have been developed ||Tl|2l|9l[l0l|27l, the generalization we introduce in this paper is, to the best of our 
knowledge, a novel approach for this problem. 
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The well known NS-Sampling Theorem IHEOlEllEllESl states that if 



and supp((7) C [— T, T] for some T > 0, then both / and g can be reconstructed from point samples of /. 
In particular, if e < 2^ then 

f{t)= 2_. /(A;e)sinc I j L^ and unif. convergence, 
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00 



g{-) = e Y^ f{ke)e^^''^- L^ convergence. 

fc=— 00 

The quantity ^, which is the largest value of e such that the theorem holds, is often refen^ed to as the 
Nyquist rate [23]. In practice, when trying to reconstruct / or g, one will most likely not be able to access 
the infinite amount of information required, namely, {/(A;e)}fegz- Moreover, even if we had access to all 
samples, we are limited by both processing power and storage to taking only a finite number. Thus, a more 
realistic scenario is that one will be given a finite number of samples {/(/>;e)}|fc|<^, for some A^ < 00, 
and seek to reconstruct / from these samples. The question is therefore: are the approximations 

fN{-) = Y. /(^e)sinc — ^ , g^{.) = eY fike)^-''^- 

k=-N \ ^ J k=-N 

optimal for / and g given the information {/(A:e)}|fc|<7v? To formalize this question consider the follow- 
ing. For A^ G N and e > 0, let 

Ojv,. = {e e C2^+' : ^ = {/(A:e)}|fc|<^, /GL2(M)nC(M)}. (1.1) 

(C(M) denotes the set of continuous functions on R). Define the mappings (with a slight abuse of notation) 



N / , , X A^ 



AiV,e,i(/)= Y. /(Msinc 1^ A7V,.,2(/) = e Y /(Me'™''- (1-2) 
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The question is, given a class of functions Q C L(IR), could there exist mappings H7v,e,i : ^N,e — ^ L 

'-'N,€,2 : ^N,e 



and Em fO '■ ^N F ^ L"^ (K) such that 



\\'^N,e,l{f) -/||l°°{R) < ||AAr,e,l(/) " /IIl°°{R) V/, / = Tg,g £ Q, 
\\^N,e,2{f) - gh'^R) < ||AAr,e,2(/) " gh^iR) V/, f = Tg,g £Q. 

As we will see later, the answer to this question may very well be yes, and the problem is therefore to find 
such mappings H7v,e,i and H7v,e,2- 

As motivation for our work, consider the following reconstruction problem. Let g be defined by 
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Figure 1: The figure shows A7v,e,2(/) for / = Fg, N = 500 and e = 0.5 (left) as well as g (right). 



This is the well-known Haar wavelet. Due to the discontinuity, there is no way one can exactly reconstruct 
this function with only finitely many function samples if one insists on using the mapping A^r ^ 2- We have 
visualized the reconstruction of g using A^r ^ 2 in Figure[T] In addition to g not being reconstructed exactly, 
the approximation KN,e,2{g) is polluted by oscillations near the discontinuities of g. Such oscillations are 
indicative of the well-known Gibbs phenomenon in recovering discontinuous signals from samples of 
their Fourier transforms ifTTl . This phenomenon is a major hurdle in many applications, including image 
and signal processing. Its resolution has, and continues to be, the subject of significant inquiry |[25l . 

However, it is tempting to think that one could construct a mapping '^N,e,2 that would yield a better 
result. Suppose for a moment that we do not know g, but we do have some extra information. In particulars 
suppose that we know that g & Q, where 
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(1.3) 



for some finite number M and where {^pk} are the Haar wavelets on the interval [0, 1]. Could we, based 
on the extra knowledge of 6, construct mappings H7v,e.i : ^N,e — ^ L^(R) and Eiv ^ •?. ■ ^iv ^ —5- L'^' 
such that 



'-'N,e,2 



SUp{||H7V,,,i(/) -/||loo( 

sup{||HAr,e,2(/) - 9\\lH 



■.gGej = Fg}< sup{||Ajv,e,i(/) - /IIl- 
: g €&,/ = Tg\ < sup{||AAr,e,2(/) - gh^i 
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Indeed, this is the case, and a consequence of our framework is that it is possible to find '^N,e,i and ^N,t,2 
such that 

sup{||Sjv,,,i(/) - /IIl^(m) : gGQJ = Tg} = 0, 
sup{||HAr,,,2(/) - 9||l2(ir) : g eej = Fg) = 0, 

provided N is sufficiently large. In other words, one gets perfect reconstruction. Moreover, the recon- 
struction is done in a completely stable way. 

The main tool for this task is a generalization of the NS-Sampling Theorem that allows reconstructions 
in arbitrary bases. Having said this, whilst the Shannon Sampling Theorem is the most standard example, 
the framework we develop addresses the more abstract problem of recovering a vector (belonging to some 
separable Hilbert space %) given a finite number of its samples with respect any Riesz basis of %. 



1.1 Organization of the Paper 

We have organized the paper as follows. In Section |2] we introduce notation and idea of finite sections of 
infinite matrices, a concept that will be crucial throughout the paper, and in Section |3] we discuss existing 
Uterature on this topic, including the work of Eldar et al |171[8l|27l. The main theorem is presented and 
proved in SectionUl where we also show the connection to the classical NS-Sampling Theorem. The error 
bounds in the generalized sampling theorem involve several important constants, which can be estimated 
numerically. We therefore devote Section [5] to discussions on how to compute crucial constants and 
functions that are useful for providing error estimates. Finally, in Section [6] we provide several examples 
to support the generalized sampling theorem and to justify our approach. 



2 Background and Notation 

Let i denote the imaginary unit. Define the Fourier transform T by 



where, for vectors a;,yGlR,x-y = a;iyi + ... + XdVd- Aside from the Hilbert space L (E ), we now 
introduce two other important Hilbert spaces: namely. 



Ql < OO 



and 
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\) = \oL = {ai,Q2,...} : ^ 
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|/3 = {.../3_i,/3o,/3i...}: J]|/32|<ooi, 
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with their obvious inner products. We will also consider abstract Hilbert spaces. In this case we will use 
the notation 'K. Note that {ej]j&i, and {ej}j^z will always denote the natural bases for /^(N) and /^(Z) 
respectively. We may also use the notation T-L for both P(N) and P{Z) (the meaning will be clear from the 
context). Throughout the paper, the symbol (g) will denote the standard tensor product on Hilbert spaces. 
The concept of infinite matrices will be quite crucial in the theory, and also finite sections of such 
matrices. We will consider infinite matrices as operators from both /^(N) to /^(Z) and /^(N) to /^(N). The 
set of bounded operators from a Hilbert space T-Li to a Hilbert space ■H2 will be denoted by B{'Hi,'H2). 
As infinite matrices are unsuitable for computations we must reduce any infinite matrix to a more tractable 
finite-dimensional object. The standard means in which to do this is via finite sections. In particular, let 
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C/gB(/2(N),/2(Z)). 



For n E N, define P„ to be the projection onto spanjei, . . . , e„} and, for odd ttt, E N, let Pm be the 



projection onto spanje m-i , . . . , C m-i }. Then PmUPn may be interpreted as 

2 2 

/ U_ m-1 1 ... U_ m-1 1 



2 ' 2 



am X n finite section of f/. Finaiiy, the spectrum of any operator T G B{'H) will be denoted by (j{T). 

3 Connection to Earlier Work 

The idea of reconstructing signals in arbitrary bases is certainly not a new idea and this topic has gone 
through extensive investigations in the last several decades. The papers by Unser and Aldroubi II1I221 have 
been very influential and these ideas have been generalized to arbitrary Hilbert spaces by Eldar ||7][8l. The 
abstract framework introduced by Eldar is very powerful because of its general nature. Our framework 
is based on similar generalizations, yet it incorporates several key distinctions, resulting in a number of 
advantages. Before introducing this framework, let us first review some of the key concepts of HI. 

Let "H be a separable Hilbert space and let f £ T-L be an element we would like to reconstruct from 
some measurements. Suppose that we are given linearly independent sampling vectors {sfcjfceN that span 
a subspace S C Ti and form a Riesz basis, and assume that we can access the sampled inner products 
Cfc = (sfe, f), k = 1,2.... Suppose also that we are given linearly independent reconstruction vectors 
{wfcjfcgN that span a subspace W C T-L and also form a Riesz basis. The task is to obtain a reconstruction 
/ G W based on the sampling data {cfcjfcgN- The natural choice, as suggested in HI, is 

f = W{S*W)-^S*f, (3.1) 

where the so-called synthesis operators S, W : P{N) — )■ ^ are defined by 

Sx = xisi + X2S2 + . . . , Wy = yiwi + y2W2 + ■ ■ ■ , 

and their adjoints S*, W* : Ti -> /^(N) are easily seen to be 

S*g = {{sug), {S2,9), . . .}, W*h = {{wi,h), {W2,h) . . .}. 

Note that S*W will be invertible if and only if 

This gives a very nice and intuitive abstract formulation of the reconstruction. However, in practice we 
will never have the luxury of being able to acquire nor process the infinite amount of samples {sk, /), 
A: = 1, 2 . . ., needed to construct /. An important question to ask is therefore: 

What if we are given only the first ?n G N samples {sk, f), k = 1, . . . , ml In this case we cannot 
use (13.11 ). and we may ask: what to do? 

Fortunately, there is a simple finite dimensional analog to the infinite dimensional ideas discussed above. 
Suppose that we are given m G N linearly independent sampling vectors {si, . . . , Sm} that span a sub- 
space Sm C T-L, and assume that we can access the sampled inner products c^ = (s^, /), fc = 1, . . . , m. 
Suppose also that we are given lineai^ly independent reconstruction vectors {wi, . . . , Wm} that span a sub- 
space Wm C T-L. The task is to construct an approximation / G Wm to / based on the samples {cfc}^^. In 



particular, we are interested in finding coefficients {dk}^^i (that are computed from tiie samples {ck}^^i) 
such that / = J2T=i (^kWk- The reconstruction suggested in |[6l is 

m 

f = Y^ dkWk = W^{S*^Wn^)-'S*^f, (3.2) 

fc=i 

where the operators Sm , Wm : C™ — ^ ^ are defined by 

SmX = XiSi + ...+XmSm, WmV = VlWi + . . . + ymWm, (3.3) 

and their adjoints S*,W* : H —?■ C™ are easily seen to be 

Sm9 = {{si,g),... ,{Sm,g)}, W^.h = {{wi,h),...,{Wm,h)}. 

From this it is clear that we can express S^Wm : C™ — )■ C" as the matrix 

{si,wi) ... (si,u;m)\ 

: : : . (3.4) 

^{Sm,Wl) ■■■ {Sm,Wm)/ 

Also, S^Wm is invertible if and only if and (lH Prop. 3]) 

WmnSi = {0}. (3.5) 

Thus, to construct / one simply solves a linear system of equations. The error can now conveniently be 
bounded from above and below by 

11/ - Pw^fW < 11/ - f\\ < —j^ ^11/ - Pw^fl 



Vjn^Jm. ) 



where PyVm i^ the projection onto Y^m, 

cosie^rnSj = mf{\\Ps„,g\\ : g G W™, H^H = 1} 

is the cosine of the angles between the subspaces Sm and Wm and Ps^ is the projection onto Sm O. Note 
that if / G Wm, then f = f exactly, a feature known as perfect recovery. Another facet of this framework 
is so-called consistency: the samples {sj, /), j = 1, . . . , m,, of the approximation / are identical to those 
of the original function / (indeed, /, as given by (13.21 ). can be equivalently defined as the unique element 
in Wm that is consistent with /). 

Returning to this issue at hand, there are now several important questions to ask: 

(i) What if Wm H Sm 7^ {0} so that SmWm. is not invertible? It is very easy to construct theoretical 
examples such that SmWm is not invertible, however (as we will see below), such situations may 
very well occur in applications. In fact, Wm H Sm = {0} is a rather strict condition. If we have 
that Wm n Sm 7^ {0} does that mean that is is impossible to construct an approximation / from the 
samples -S*^/? 

(ii) What if IK-SmVFm)"^!! is large? The stability of the method must clearly depend on the quantity 
||(S'*„Wm)~^||- Thus, even if (5j^Wm)~^ exists, one may not be able to use the method in practice 
as there will likely be increased sensitivity to both round-off eiTor and noise. 
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Figure 2: This figure shows logio ll(5'*m^m)~^|| as a function of m and e for m = 1, 2, . . . , 100. The 
left plot corresponds to e = 1, whereas the right plot corresponds to e = 7/8 (circles), e = 1/2 (crosses) 
and e = 1/8 (diamonds). 



Our framework is specifically designed to tackle these issues. But before we present our idea, let us 
consider some examples where the issues in (i) and (ii) will be present. 

Example 3.1. As for (i), the simplest example is to let Ti = /^(Z) and {ej}j£z be the natural basis (cj 
is the infinite sequence with 1 in its j-th coordinate and zeros elsewhere). For m € N, let the sampling 
vectors {sfc}^_,^ and the reconstruction vectors {t(;fc}^_^ be defined by Sk = e^ and Wk = eu+i- 
Then, clearly, Wm n 5^ = span{em+i}- 

Example 3.2. For an example of more practical interest, consider the following: For < e < 1 let 

T-L = -L^([0, 1/e]), and, for odd m G N, define the sampling vectors 



l*e,KJfc=-(m-l)/2' 



Se,fc 



e X[0,l/e]' 



(this is exactly the type of measurement vector that will be used if one models Magnetic Resonance 
Imaging) and let the reconstruction vectors {it'/d^i denote the m first Haar wavelets on [0, 1] (including 
the constant function, wi = X[o,i])- Let Se^m and Wm be as in ( 13.31 ). according to the sampling and 
reconstruction vectors just defined. A plot of || (5'*^M^m)~^ || as a function of m and e is given in Figure 
|2] As we observe, for e = 1 only certain values of m yield stable reconstruction, whereas for the other 
values of e the quantity ||(<S'*„^Wm)~^|| grows exponentially with m, making the problem severely ill- 
conditioned. Further computations suggest that ||(<S'*^VFm)~^|| increases exponentially with m not just 
for these values of e, but for all < e < 1. 

Example 3.3. Another example can be made by replacing the Haar wavelet basis with the basis consisting 
of Legendre polynomials (orthogonal polynomials on [—1, 1] with respect to the Euclidean inner product). 
In Figure [3]we plot the quantity ||(5'*„Wm)~^||. Unlike the previous example, this quantity grows 
exponentially and monotonically in m. Whilst this not only makes the method highly susceptible to round- 
off en^or and noise, it can also prevent convergence of the approximation / (as m — )• oo). In essence, for 
convergence to occur, the en^or ||/ — -P>v„j/|| must decay more rapidly than the quantity IKS*^^^^)"^!! 
grows. Whenever this is not the case, convergence is not assured. To illustrate this shortcoming, in Figure 
|3]we also plot the error ||/ — /||, where /(x) = x+i6z^ • "^^^ complex singularity at x = ibii hmits the 
convergence rate of ||/ — Pw„J\\ sufficiently so that / does not converge to /. Note that this effect is 
well documented as occurring in a related reconstruction problem, where a function defined on [— 1 , 1] is 
interpolated at m equidistant pointwise samples by a polynomial of degree m — 1. This is the famous 
Runge phenomenon. The problem considered above (reconstruction from m Fourier samples) can be 
viewed as a continuous analogue of this phenomenon. 
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Figure 3: The left figure shows log^o || {S*„^Wm) ^ || as a function of rra for ?ti = 2, 4, . . . , 50 and e = 
1, |, i, I (squares, circles, crosses and diamonds respectively). The right figure shows log^Q ||/ — Pw^fW 
(squares) and logio ||/ - f\\ (circles) for m = 2, 4, 6, ... , 100, where f{x) = j^jq^- 

Actually, the phenomenon illustrated in Examples l3.2l and [3.3l is not hard to explain if one looks at the 
problem from an operator-theoretical point of view. This is the topic of the next section. 



3.1 Connections to the Finite Section Method 

To illustrate the idea, let {sfcjfceN and {uifcjfcgN be two sequences of Unearly independent elements in a 
Hilbert space T-L. Define the infinite matrix U by 



U 



I'uii Ui2 Ui3 . . .^ 
U2l U22 U23 ■ ■ 
U3l U32 U33 . . 



\ 



Uij = {Si,Wj). 



(3.6) 



■/ 



Thus, by (13.41 ) the operator S^Wm is simply the m by m finite section of U. In particular 



S*mW„^ 



P UP 



The finite section method has been studied extensively for the last decades |l3l[l2l[T3l|2Tl. It is well known 
that even if U is invertible then Pm.UPm\p^i^(m may never be invertible for any m. In fact one must have 
rather strict conditions on U for PmUPmlp^i^m) to be invertible with uniformly bounded inverse (such 
as positive self-adjointness, for example II2TI ). In addition, even if U : P{N) — )■ P(N) is invertible and 
Pm.UPm\prni'^m) is invertible for all tti G N, it may be the case that, if 



X = U ^y, 



X,y £ /^(N), Xr,^ = {PmUPrn\p^l2(N)) ^ PmV 



then 






m 
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Suppose that {sfcjfcgN and {t(;fc}fceN are two Riesz bases for closed subspaces S and lA^ of a separable 
Hilbert space li. Define the operators 5, W -.l^i^) ^li by 



Sx = xisi + X2S2 + — 



Wy = yiwi + y2W2 + .... 



(3.7) 



Suppose now that {S*W)~^ exists. For tti G N, let the spaces 5^, Wm and operators Sm, Wm. '■ C"^ — )• Ti 
be defined as in Section |3] according to the vectors {sk}^^i and {wk}^^i respectively. The following 
scenarios may now arise: 



(i) WnS^ ={0},yet 

WmnS^^iO}, VmGN. 

(ii) IKS*!^)"^!! < c» and the inverse {S^Wm)~^ exists for all m G N, but 

\\{S*^Wm)~^ \\ ^ oo, m^oo. 

(iii) {S*^Wm)~^ exists for all m G N, however 



Wm{S*^Wmr'SU ^ /, m ^ oo, 



for some / G W. 



Thus, in order for us to have a completely general sampling theorem we must try to extend the frame- 
work described in this section in order to overcome the obstacles listed above. 



4 The New Approach 

4.1 The Idea 

One would like to have a completely general sampling theory that can be described as follows: 

(i) We have a signal / G "H and a Riez basis {wfcjfegM that spans some closed subspace W CH, and 



f = Y^hwk, /3fe G C. 



fc=i 

So / G W (we may also typically have some information on the decay rate of the /3fcS, however, 
this is not crucial for our theory). 

(ii) We have sampling vectors {sfcjfcgN that form a Riez basis for a closed subspace S C T-L, (note that 
we may not have the luxury of choosing such sampling vectors as they may be specified by some 
particular model, as is the case in MRI) and we can access the sampling values {(s^, /)}fceN- 

Goal: reconstruct the best possible approximation / G W based on the finite subset {{sk, /)}fcLi of the 
sampling information {(s^, /)}fc6N- 

We could have chosen m vectors {wi, . . . , Wm} and defined the operators Sm and Wm as in (13.31 ) 
(from {wi, . . . , Wrn} and {si, . . . , Sm}) and let / be defined by ( 13.21 ). However, this may be impossible 
as S'^Wm may not be invertible (or the inverse may have a very large norm), as discussed in Examples 
I3.2l and l331 The question is then: what to do? 

To deal with these issues we will launch an abstract sampling theorem that extends the ideas dis- 
cussed above. To do so, we first notice that, since {sj} and {wj} are Riesz bases, there exist constants 
A,B,C,D > such that 



A}]\a,\'< 
fceN 


y ]akWk 
fceN 
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y]akSk 

fceN 
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fceN 

<DY,Wk\' 

k&i 



(4.1) 



V{Qi,a2,...}G/2(N). 



Now let U be defined as in (13.61) . Instead of dealing with PmUPm\p^i2{n) = S^Wm we propose to 
choose n G N and compute the solution {/3i , . . . , /3„} of the following equation: 



A 



//3i\ 



P U*P 



{S2j) 



\{Sm,f)J 



A = PnU*P^UPJp^^ 



n, 



(4.2) 



provided a solution exists (later we will provide estimates on the size of n, m for ( 14.21 ) to have a unique 
solution). Finally we let 

n 

f = Y,hwk. (4.3) 

fc=i 

Note that, for n = m this is equivalent to (13.21) . and thus we have simply extended the framework discussed 
in Section [3l However, for m > n this is no longer the case. As we later establish, allowing m to range 
independently of n is the key to the advantage possessed by this framework. 

Before doing so, however, we first mention that the framework proposed above differs from that 
discussed previously in that it is inconsistent. Unlike ( 13.21 ). the samples {sj, f) do not coincide with those 
of the function /. Yet, as we shall now see, by dropping the requirement of consistency, we obtain a 
reconstruction which circumvents the aforementioned issues associated with (13.21). 



4.2 The Abstract Sampling Theorem 

The task is now to analyze the model in (14.21) by both estabhshing existence of / and providing error 
bounds for 11/ — /||. We have 

Theorem 4.1. Let Ti be a separable Hilbert space and S,W C Ti be closed subspaces such that yV n 
S = {0}. Suppose that {sk}k&l <^nd {wk}keN <^re Riesz bases for S and W respectively with constants 
A,B,C,D > 0. Suppose that 



f = Y.PkWk, /3 = {/3i,/32,...,}G/2(pj)_ 



(4.4) 



Let n € N. Then there is an M € N (in particular M = ininjA; : ^ a{PnU*PkUPn\p^-}{)}) such that, 
for all m > M, the solution {/3i , . . . , /3„} to i \4.2i is unique. Also, if f is as in ( 14. J D . then 



where 



Kr. 



||/-/||w<VS(l + if„,„)||P„^/3||p(N), 



{PnU*P^UPn\p„n)-^PnU*PmUP,i 



(4.5) 
(4.6) 



The theorem has an immediate corollary that is useful for estimating the error. We have 
Corollary 4.2. With the same assumptions as in Theorem \4.1\ and fixed n G N, 

1 



iPnU*Pn^UPn\p„n)~^ 



{PnU*UPn\p„nr^\ < \\iU*Ur^ <—, m^oo. 



(4.7) 



In addition, ifU is an isometry (in particular, when {wk}keNi {sk}keN cire orthonormal) then it follows 
that 



Kr, 



0, 



m — 7- 00. 
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Proof of Theorem \4J\ Let U be as in as in (13.61) . Then (14.41) yields the following infinite system of equa- 
tions: 

f{siJ)\ All Ui2 Ui3 ...\ //3i^ 

(S2,/) U21 U22 U23 ••• /32 

(■S3,/) ^ U3I U32 U33 ... fSs 



\ ■■ J {■■ 



(4.8) 



•••/ \-J 



Note that U must be a bounded operator. Indeed, let S and W be as in (13.71 ). Since 

(5*VFej,ei) = (si,u;j), i, j G N, 

it follows that U = S*W. However, from (14.11 ) we find that both W and S are bounded as mappings from 
^^(N) onto W and S respectively, with ||VF|| < Vb, \\S\\ < VD, thus yielding our claim. Note also that, 
by the assumption that W n S^ = {0}, (14.81 ) has a unique solution. Indeed, since W n 5-*- = {0} it 
follows that inf|j^||=i ||S'*14/^x|| / 0, so [/ must be injective. 
Now let r]f = {(si, /), (si, /), . . .}, then (1481 ) gives us that 

(4.9) 



PnU*PmVf = PnU*PmU P„ + P„^ /3 



Suppose for a moment that we can show that there exists an M > such that PnU* PmU Pn\p„H is 
invertible for all m > M. Hence, we may appeal to ( 14.91 ), whence 



(4.10) 



{PnU*Pm.UPn\p„nr'PnU*P„,Vf = PnP + {PnU* P^U Pn\p^nr' PnU* P„,.UP^ 13, 

and therefore, by (l4!9l ) and (1411 ). 

<^\\{PnU*PmUPn\p,^nr'PnU*Pmr]f-p\L(„. 



f -^ hwk 



k=l 



n 



{P^ - iPnU*PmUPn\p„n)~^PnU*PmUP;t)P 



«2(N) 



<^{l + Kn,m) Pn^ 



Z2(N) 



where 



K„ 



:PnU*P„,UPn\pM)~^PnU*PmUP;i 



Thus, (14.51 ) is established, provided we can show the following claim: 

Claim: There exists an M > such that PnU*PmUPn\p„n is invertible for all m > M. Moreover, 



\{PnU*PmUPn\p„nr^\ 



{PnU*UPn\p^n)~^\\< \\iU*U)'H, 



m — )• oo. 



To prove the claim, we first need to show that PnU*UPn\p^i2m\ is invertible for all ?i G N. To 
see this, let G : B{P{'N)) — )• C denote the numerical range. Note that U*U is self-adjoint and invert- 
ible. The latter implies that there is a neighborhood a; around zero such that a{U*U) H u = and the 
former implies that the numerical range @{U*U) n w = 0. Now the spectrum a{PnU*UPn\p^i2m^) C 
e{PnU*UPn\p„i2(^n)) C e{U*U). Thus, 



a{PnU*UPn\p^i2^n))ri 



UJ 



VnGN, 
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thus, PnU*UPn\p^p(n) is always invertible. Now, make the following two observations 



PnU*UPn 






(4.11) 






(^ne 



J^ 



where the last series converges at least strongly (it converges in norm, but that is a part of the proof). The 
first is obvious. The second observation follows from the fact that Pm.U — )■ U strongly as m — )• oo. Note 
that 

However, U*PnU must be trace class since ran(P„) is finite-dimensional. Thus, by (14.111 ) we find that 



\Pr,U*PrnUPr, " Pr,U*U PA < 



\{PnCj) (^ {PnCj)\\ 



oo 

E 

j=m+l 

oo 

< Yl {UPnU*ej,ej, 

j=m+l 



0, 



m 



oo. 



Hence, the claim follows (the fact that II (P„C'^*C'^P„ I P„-^) ^|| < ||(t/*C/) ^ 1 1 is clear from the observation 
that U* U is self-adjoint) and we are done. D 

Proof of Corollary 14. 21 Note that the claim in the proof of Theorem 14. 1 1 yields the first part of (14.71 ). and 
the second part follows from the fact that U = S*W (where 5, W are also defined in the proof of Theorem 
14.11 ) and (14.11 ). Thus, we are now left with the task of showing that Kn,m — )■ as m — ;• oo when U is an 
isometi"y. Note that the assertion will follow, by (14.61) . if we can show that 



0, 



m 



oo. 



P UP UP 
However, this is straightforward, since a simple calculation yields 

< \\U\\{\\PnU*PmUPn " PnU*UPn 



P U*P UP 



a/2 



(4.12) 
D 



Remark 4.3. Note that the trained eye of an operator theorist will immediately spot that the claim in the 
proof of Theorem [4.11 and Corollary l4.2l follows (with an easy reference to known convergence properties 
of finite rank operators in the strong operator topology) without the computations done in our exposition. 
However, we feel that the exposition illustrates ways of estimating bounds for 



{PnU*P^UPn\p„ny^ 



Pr,U*PmUP' 



which are crucial in order to obtain a bound for Kn^m- This is demonstrated in Section|5l 

Remark 4.4. Note that 5*1^ (and hence also U) is invertible if and only ifH = W ® S-^, which is 
equivalent to lA^ n S-^ = {0} and W'^ n 5 = {0}. This requirement is quite strong as we may very 
well have that W ^ T-L and 8 = 1-1 (e.g. Example 3.2 when e < 1). In this case we obviously have that 
W"*" n 5 / {0}. However, as we saw in Theorem 14. 1[ as long as we have / G >V we only need injectivity 
of C/, which is guaranteed when W r\S-^ = {0}. 
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If one wants to write our framework in the language used in Section |3l it is easy to see that our 
reconstruction can be written as 

/ = Wn{W;:SmS*,,Wn)~^W:^SmS*^f, (4.13) 

where the operators S^ : C" — )■ 7i and Wn : C" — )■ H are defined as in ( I3.3I ). and Sm and Wn 
corresponds to the spaces 

Sm = span{si,...,Sm}, Wn = span{'u;i,...,u;„}, (4.14) 

where {tiifcjfceN and {sfcjfceN are as in Theorem 14.11 In particular, we get the following corollary: 

Corollary 4.5. Let T-Lbe a separable Hilbert space and 5, W C T-Lbe closed subspaces such that W n 
S^ = {0}. Suppose that {sfcjfcgN <^nd {ii^ytlyteN ^^^ Riesz. bases for S and W respectively. Then, for each 
n G N there is an M gN such that, for all m > M, the mapping W^SmS^Wn : C^ — t- C" is invertible 
(with S„i and Wn defined as above). Moreover, if f is as in d^.iJD . then 



P^J ^, < 11/ - fWn < (1 + Kn,m) P^J 



H 



n 



where Py\!^ is the orthogonal projection onto Wn, and 



Kr, 



Wn{W:,SmS*n,Wn)-^W*Sn,S*„,P<^^ 



Moreover, when {sfc} and {wk} are orthonormal bases, then, for fixed n, Cn,m — ?• a^ m — )■ oo. 

Proof. The fact that WnSmS^Wn : C" — )• C" is invertible for large m follows from the the observation 
that W nS-^ = {0} and the proof of Theorem HU by noting that S^Wn = PmUPn, where U is as in 
Theorem 14.11 Now observe that 



W::SmS*n,f = W^Sn^SUPwJ + PwJ) 

= W:SmS*n,Wn{W:Wnr'W:j + W*Sn,S*n,P^J. 



(4.15) 



Note also that W*Wn:C 
(14.151 ) yields 



> C"- is clearly invertible, since {wk}^^i are hnearly independent. Now 

Wn{W*SmS*n,Wn)-'W:SmS*„,f = P^J + Wn{W:SmS*n,Wn)-'W:SmS*n,P^J. 



Thus, 



\\f-f\\n< 



P^„-Wn{W:SmS*n,Wnr'W:SmS*n,P^Jv,, P^J 

H 



n 



which gives the first part of the corollary. The second part follows from similar reasoning as in the proof 
of Corollary SH □ 

Remark 4.6. The framework explained in Section [3] is equivalent to using the finite section method, and 
this may work for certain bases, however, not in general (as Example l3.2l shows). Computing with infinite 
matrices can be a challenge since the qualities of any finite section may be very different from the original 
infinite matrix. The use of uneven sections (as we do in this paper) of infinite matrices seems to be the 
best way to combat these problems. The reader may consult ifTTl [141 [151 for other examples of uneven 
section techniques. 
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When compared to the method of Eldar et al, the framework presented here has a number of important 
advantages: 



(i) It allows reconstructions in arbitrary bases and does not need extra assumptions as in (13.51) . 

(ii) The conditions on m (as a function of n) for PnU* PmU Pn\p„n to be invertible (such that we have 
a unique solution) can be numerically computed. Moreover, bounds on the constant Kn,m can also 
be computed efficiently. This is the topic in Section [5] 

(iii) It is numerically stable: the matrix A = PnU*PmUPn\p„n has bounded inverse (Corollary 14.21 ) for 
all n and m sufficiently lai^ge. 

(iv) The approximation / is quasi-optimal (in n). It converges at the same rate as the tail ||P^/3||;2™), 
in contrast to ( 13.21) which converges more slowly whenever the parameter — t^ r grows with 

n = m. 

As mentioned, this method is inconsistent. However, since {sj} is a Riesz basis, we deduce that 

m 

Y.\{s„f-f)\'<c\\f-ff, 

for some constant c > 0. Hence, the departure from consistency (i.e. the left-hand side) is bounded by a 
constant multiple of the approximation error, and thus can also be bounded by ||P^/3||;2m). 

4.3 The Generalized (Nyquist- Shannon) SampUng Theorem 

In this section, we apply the abstract sampling theorem (Theorem 14.11 ) to the classical sampling problem 
of recovering a function from samples of its Fourier transform. As we shall see, when considered in this 
way, the corresponding theorem, which we call the generalized (Nyquist-Shannon) Sampling Theorem, 
extends the classical Shannon theorem (which corresponds to a special case) by allow reconstructions in 
arbitrary bases. 

Proposition 4.7. Let J- denote the Fourier transform on Lp'iMr). Suppose that {<^j}jgN is a Riesz basis 
with constants A, B (as in ( 14.71 )) for a subspace W C L'^{W^) such that there exists a T > with 
supp{ipj) C {—T,T]'^ for all j G N. For e > 0, let p : N ^ [eLY be a bijection. Define the infinite 
matrix 

(Uii Ui2 Ui3 . . .\ 



U 



U21 U22 U23 
U31 U32 U33 



Uij = {Tipj){p{i)). 



V : : : ■•/ 

Then, for e < 1^, we have that U : /^(N) — >■ i^(N) is bounded and invertible on its range with \\U\\ < 
Ve^ and \\{U*U)~^\\ < e'^A^^ . Moreover, if {^j}j^n is an orthonormal set, then e'^'^C/ is an 
isometry. 

Theorem 4.8. (The Generalized Sampling Theorem) With the same setup as in Proposition \4.7\ set 

00 

14 



and let P^ denote the projection onto spanjei , . . . , e^}- Then, for every n G N there is an M £ 'N such 
that, for all m > M, the solution to 



A 



//3i\ 

V/3n/ 



P JJ*P 



n'-^ J m. 



//(P(l))\ 

/(P(2)) 
\fipim))) 



A = PnU*PmUPn\p„n, 



is unique. Also, if 

then 
and 



9 = ^^jVj, 






\9 - 9\\lHR^) < ^(1 + Kn,rn)\\PnP\\pm^ ^ = i/^l'/^a, • • •}, 



(4.16) 
(4.17) 



where Kn^m i^ given by Ii4.6i and satisfies ( 14. 71) . Moreover, when {ipj}j£^ is an orthononnal set, we have 

for fixed n. 

Proof of Proposition 14. 71 Note that 



ra — )• oo, 



Ui^= / <^j(x)e-2"'''(^)-^dx 



'T^TY 



^j (x) e-'^'^'^P^^-^ dx. 



Since p : N — ;■ (eZ)^ is a bijection, it follows that the functions {x i— )■ e'^/^e ^'^^^(^^'^jjgN form an 
orthonormal basis for L2([-(2e)-^(2e)-l]'^) D L2([-r,T]'='). Let 



(•,•) = (•,•) 



L2([-(2.)-i,(2.)-i]d), 



denote a new inner product on L^([— (2e)^^, (2e)~^]'^). Thus, we are now in the setting of Theorem 14. 11 
and Corollary 14.21 with C = D = e'^. It follows by Theorem 14. 1 1 and Corollary 14.21 that U is bounded 
and invertible on its range with \\U\\ < ve^ and ||(C/*C/)~^|| < e'^A"^. Also, e'^-^'^U is an isometry 
whenever A = i3 = 1, in particular- when {c^fcjfceN is an orthonormal set. D 



Proof of Theorem \4~8\ Note that (14.161) now automatically follows from Theorem 14. II To get (14.17b we 
simply observe that, by the definition of the Fourier transform and using the Cauchy-Schwarz inequality. 



sup 






< 



9iy) -^k'Pjiy) 



dy 



< (2T) 



d/2 



9-Yl ~^i^i 






L2 



where the last inequality follows from the already established (14.161 ). Hence we are done with the first 
part of the theorem. To see that Kn^m — 5- as m — ;• oo when {ipj^j^-^ is an orthonormal set, we observe 
that orthonormality yields A = B = 1 and hence (since we already have established the values of C and 
D) e'^/'^U must be an isometry. The convergence to zero now follows from Theorem 14. II D 
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Note that the bijection p : N — > {eL^ i^ oi^ly important when d > 1 to obtain an operator U : /^(N) — )• 
P(N). However, when d = I, there is nothing preventing us from avoiding p and forming an operator 
U : P{N) — )■ /^(Z) instead. The idea follows below. Let T denote the Fourier transform on L^(]R), and 
let / = Tg for some g G L^(M). Suppose that {^j}jeN is a Riesz basis for a closed subspace in L^(M) 
with constants A,B>0, such that there is a T > with siipp{(pj) C [-T, T] for all j G N. For e > 0, 
let 



U 



( : : : 

"-1,1 lf-1,2 ""-1,3 

"0,1 "0,2 "0,3 

^1,1 ^^1,2 "1,3 

V : : ; 



Ui^j = {Fipj){ie). 



(4.18) 



Thus, as argued in the proof of Theorem WM U G B{P{'N),P{Z)), provided e < ^. Next, let P„ G 
B{P{N)) and, for odd m, Pm G B{P{Z)) be the projections onto 

spanjei, . . . ,6^}, span{e_m-i , . . . ,em-i } 

2 2 

respectively. Define {/5i, . ..,/?„} by (this is understood to be for sufficiently large m) 



A 



h 



P U*P 



n^ ^ m 



/(O) 



V /(^) / 



) A — PnU* PmJJ Pn\p,M- 



(4.19) 



By exactly the same arguments as in the proof of Theorem 14.81 it follows that, if g = YlTLil^j^j^ 

9 = Ej=i ^jfj' f = J='g and f = Yl]=i i^j^Vj, then 



\\g - 5||l2(k) < /B(1 + ir„,„)||P„^/3||p(j,), /3 = {/3i,/32, • • •}, 

11/ - /IIloo(k) < V2fVB{l + Kn,m)\\P;t ^\\piN), 

where Kn m is as in (I4.6I ). 



(4.20) 



Remark 4.9. Note that (as the proof of the next corollary will show) the classical NS-Sampling Theorem 
is just a special case of Theorem 14. 8 1 

Corollary 4.10. Suppose that f = Fg and supp(g) C [— T, T]. Then, for {) < e < lyp we have that 



ff(-) = e E f{ke)e 



L convergence. 



k=—oo 



fit) = y f{ke)smc I ) L and unif. convergence. 



k=—oo 
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Proof. Define the basis {(/JjjjeN for L2([-(2e)-\ (2e)~i]) by 



<- 2e'2eJ L 2£'2£J 

t 2e'2eJ 



(^6(x) = V^e^™'^^X[_^,^](2;) etc. 
Letting U = {u^ /Ifcgz^^gN, where n^ / = {Tipi){ke), an easy computation shows that 



/ 



U 















i 








1 

7^ 








1 

















1 




















1 






v^ 



V : 



By choosing m = n in (14.191 ). we find that /3i = ^/ef{0), f32 = \/e/(e), ^3 = ^/ef{-e), etc and that 
Kn,m = in (14.201 ). The corollary then follows from (14.20b . D 

Remark 4.11. Returning to the general case, recall the definition of QN,e from (II. lb . the mappings AN,e,i, 
Km fO from (11.21) and G from (11.31). Define Hat^ 1 : ilAr^ — ^ L^(M) and Hat^ ;> : ilAr^ — ^ L^' 



-AT, 



-Af,e,l 
n 



-Af,e,2 : i'Af,( 
n 



by 



where /3 = {f3i, . . . , /3.„} is the solution to ( 14.191 ) with N = m. Then, for n > M (recall M from the 
definition of G ( 11.31 )). and 



m 



7n(7) = min{fc G N : ||(P„^*Pfct/P„|p„w)-i|| < £7}, 7 > 1, 



it follows that 



W^N.eM) -/||l°o(k) =0 < ||AAr,e,l(/) - f\\L^[K)^ f J = J^Q, 9 G Q, 
pAr,e,2(/) " 5||l2(R) = < || AAr,e,2(/) " 5||l2(R)V/, f = Fg,g ^Q. 

Hence, under the aforementioned assumptions on m and n, both / and g are recovered exactly by this 
method, provided (7 G Q. Moreover, the reconstruction is done in a stable manner, where the stability 
depends on the parameter 7. 



To complete this section, let us sum up several of the key features of Theorem 14.81 First, whenever 
m is sufficiently lai^ge, the error incun^ed by g is directly related to the properties of g with respect to the 
reconstruction basis. In particular, as previously noted, g is reconstructed exactly under certain conditions. 
Second, for fixed n, by increasing m we can get arbitrarily close to the best approximation to g in the 
reconstruction basis whenever the reconstruction vectors are orthonormal. Thus, provided an appropriate 
basis is known, this procedure allows for near-optimal recovery. The main question that remains, however, 
is how to guarantee that the conditions of Theorem 14. 8 1 are satisfied. This is the topic of the next section. 
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5 Norm Bounds 



5.1 Determining m 

Recall that the constant Kn^m in the error bound in Theorem l4. ll Crecall also U from the same theorem) is 
given by 



K„ 



{PnU*P„^UPn\p„nr^PnU*FmUP,f 



It is therefore of utmost importance to estimate Kn^m- This can be done numerically. Note that we 
already have established bounds on \\U\\ depending on the Riesz constants in (14.11 ) and since we obvi- 
ously have that Kn^^ < ll(-fnf^*-Pm^-fn|p„'H)~^lll|f^lP) we only require an estimate for the quantity 

\\{PnU*P^UPn\pMr\ 

Recall also from Theorem 14. 11 that, if U is an isometry up to a constant, then i^„^m — >• as m ^ oo. 
In the rest of this section we will assume that U has this quality. In this case we are interested in the 
following problem: given n G N, 6* G M+, what is the smallest ?n G N such that Kn^m < ^? More 
formally, we wish to estimate the function $ : ZY(/^(N)) x N x M+ -> N, 



$(C/,n,0) = min{mGN: ||(P„C/*P^C/P„|p„^)^ip„C/*P„C/P„^|| < ^j , (5.1) 



where 



U{f{n)) = {[/ G fi(/2(N)) :[/*[/ = c/, c G M+} . 
Note that <I> is well defined for all 6 G M+, since we have established that Kn m 



as 771 — )■ oo. 



5.2 Computing Upper and Lower Bounds on Kn,m 

The fact that UP^ has infinite rank makes the computation of Kn^m a challenge. However, we may 
compute approximations from above and below. For M G N, define 



K„ 



{PnU*PmUPn\p„n)~'PnU*PmUP;tPM 
= \\{PnU*PmUPn\p,^Hr^PnU*Pj\ . 



Then, for L > M 



sup 
< sup 



< sup 

?6-H,||f||=l 



{PnU*P^UPn\pMr^PnU*PmUP^PMi 
iPnU*PmUPn\p„n)-^PnU*PraUP;tPL^ 
PnU*PmUPn\p„nr^PnU*PmUP,iC 



K„ 



Clearly, Kn^m < ||^||-f^n,m and, since PmS, — )• ^ as M — )• oo for all ^ G Ti, and by the reasoning above, 
it follows that 

Kn,m,M < Kn,m < ||t^||-^n,m, Kn,m,M Z' ^n,m, M —?■ OO. 

Note that 

{PnU*PmUPn\pM)~^PnU*PrnUP,iPM : PmU ^ PnU 

has finite rank. Therefore we may easily compute Kn,m,M- In Figured we have computed Kn,rn,M for 
different values of n, m, M. Note the rapid convergence in both examples. 
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Figure 4: The figure sliows Kn,m,M for n = 75, m = 350 and M = n + 1, . . . , 6000 (left) and Kn,m,M 
for n = 100, m = 400 and M = n + 1, . . . , 6000 (rigiit) for tlie Haai" wavelets on [0, 1]. 

5.3 Wavelet bases 

Whilst in the general case ^{U, n, 9) must be computed numerically, in certain cases we are able to derive 
explicit analytical bounds for this quantity. As an example, we now describe how to obtain bounds for 
bases consisting of compactly supported wavelets. Wavelets and their various generalizations present an 
extremely efficient means in which to represent functions (or, more commonly, signals) ||4l|5l|22l. Given 
their long list of applications, the development of wavelet-based reconstruction methods is naturally a 
topic of utmost importance. 

Let us review the basic wavelet approach on how to create orthonormal subsets {ipk}keN C -L^(]R) 
with the property that L'^{[0, a]) C cl(span{(/9fc}fcgi^) for some a > 0. Suppose that we are given a mother 
wavelet V' and a scaling function such that supp('0) = supp((/>) = [0, a] for some a > 0. The most 
obvious approach is to consider the following collection of functions: 

^a = {(l)k,i^j,k ■.jeZ+,k£Z, supp(0fc)" n [O, a] / 0, supp{ijj^k)° n [O, a] / 0}, 

where 

ct)k = </>(• - k), ^j, fc = 2iV(2^' • -A;). 

(The notation K° denotes the interior of a set J^ C M.) Then we will have that 

L2([0,a]) C cl(span{(^ : ip G fla}) C L'^[-T,T], 

where T > is such that [—T,T] contains the support of all functions in Q,a- However, the inclusions 
may be proper (but not always, as is the case with the Haar wavelet.) It is easy to see that 

a + k k 

i^j,k ^ ^a ^^ -^— < 0, a< —, 

(f'k ^ ^a *^==^ a + k <0, a < k. 
Hence we get that 

fia = {<t>k : |A:| = 0, . . . , [a] - 1} U {Vj.fc : j e Z+,k e Z,-\a] + 1 < k < 2^ \a] - 1}, 

and we will order Q,a as follows: 

{(j),(j)l, . . . ,(j)^a\-l,<P-l, ■ ■ ■ ,4'-\a\+l,1pO,0,1pO,l, ■ ■ ■ ,'(/'0,ral-l,V'0-l,---,^0-ral+l'^l,0,---}- (5-2) 

We will in this section be concerned with compactly supported wavelets and scaling functions satisfying 

\J'(t>{w)\<j^, |-^^HI<^, a;GM\{0}, (5.3) 
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for some 



C>0, 



p G N. 



Before we state and prove bounds on ^{U,n,6) in this setting, let us for convenience recall the result 
from the proof of Theorem 14.11 In particular, we have that 



\PnU*PmUPn-PnU*UPn\\< J^ {UPnU*e,,e 

j=m+l 



Oh 



m — 7- oo. 



(5.4) 



Theorem 5.1. Suppose that {ipi}i^f^ is a collection of functions as in f l5.2D such that supp(c/3i) C [— T, T] 
for all I ^ fi and some T > 0. Let U be defined as in Proposition \4.7\ with < e < ^^ and let the 
bijection /9 : N -^ eZ defined by p(l) = 0, p(2) = e, p{3) = -e, p(4) = 2€, . . .. For 9 > 0,n e N define 
^{U, n, 6) as in Ii5.1\) . Then, ifcj), ip satisfy d5.JD , we have that 



;^,TT m/ /4ei-^''M'^^\^^ /-, /4Pn2P-l\\5^ r^ ( ^^ 



/W 



I- 



4P_1 



n — )• oo, 



where f{e) = (Vl + 4^2 _ i)2/(4^2)_ 

Proof. To estimate ^{U,n,6) we will determine bounds on 



p Tr*p Tjp^ 



< 



'}^ 



Note that if r < 1 and ||P„[/*P„[/P„ - PnU*UPn\\ < r, then \\{PnU*PmUPn\p„ny^\\ < e/(l - er) 
(recall that [/*[/ = e'^I and that e < 1). Also, recall (14.121 ). so that 



\\{PnU*P„^UPn\p^Hr' 

when r and m are chosen such that 



P U*P UP 



< 



<e. 



1 — er 
(recall that \\U\\ = l/y/e). In particular, it follows that 



PnU*PnUPn " PnU*UPn\\ < r, 



■^{U,n,e)<mm{m: \\PnU* P^U Pn - PnU* U Pn\\ < e'^ {^/l + Aff^ - if / {Ae"^)} . (5.5) 

To get bounds on ^(f/, n, 9) we will proceed as follows. Since (/>, -0 have compact support, it follows that 
F4>, Tip are bounded. Moreover, by assumption, we have that 



\:fcP{w)\ < 



c 



l-^V'HI < 



c 



w 



u gM\{0}. 



And hence, since 
we get that 



T^jj^kiw) = e~'^''''^'"2 2T^{2^iw), 



|J^Vi,fc(^)l < 2^ 



C 

\2-3w\P' 



uj G 



(5.6) 
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Note that, by the definition of U, it follows that 



oo n 



j=m+l s=m+l t=l 

And also, by ( I5.6I ) and (15.21 ) we have, for s > 0, 

n Llog2WJ 2^N-1 

t=l i=0 k=-\a]+l 



2[a]C'^ U°S2(")J 2JN-1 . C^ / C^ Lios2('^)J ^2 



2ralC2 / 4Pn^P-l 

< , '. ' o,., 1 + 



|p(s)|2p V 4P-1 

thus we get that 

4Pn2p - 1 



3=m+l t=l ^ ^ s=m+l '^^ '' 

II V 4P - 1 / -^ s2p - m2P"i V 4P - 1 

^ -' s=m+l ^ 



(5.7) 



Therefore, by using (15.41 ) we have just proved that 



\pu*pup -p^-^p||< '^"''^"^^7 i + ''"''-^ 

m^p-^ \ 4P — 1 



and by plugging this bound into (15.51) we obtain 

which obviously yields the asserted bound on ^{U, n, 6). D 

The theorem has an obvious corollary for smooth compactly supported wavelets. 

Corollary 5.2. Suppose that we have the same setup as in Theorem 15.71 and suppose also that </>, ^ G 

CP{M.) for some p G N. Then 

^{U,n,e) = o(n'^\ , n -^ 00. 



5.4 A Pleasant Surprise 

Note that if if) is the Haar wavelet and = X[o,i] we have that 

\:f<p{w)\ < ^, i^vMI <A, w G 

\w\ Iw] 
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Figure 5: The figure shows sections of the graphs of "^{U, •, 1) (left) and "^{U, •, 2) (right) together with 
the functions (in black) x i— )■ 4.9x (left) and x i— )■ 4.55x. In this case U is formed by using the Haar 
wavelets on [0, 1]. 



Thus, if we used the Haar wavelets on [0, 1] as in Theorem 15.11 and used the technique in the proof of 
Theorem 15. II we would get that 



min{7n : ||P„C/*/'^C/P„-P„t/*C/P„|| = e-\VT+49^ -1)^ /{AO^)} = O (n^) , 



n 



oo. 



(5.8) 



However, it is tempting to check numerically whether this bound is sharp. Let us denote the quantity in 
(15.81 ) by ^{U, n, 9), and observe that this can easily be computed numerically. Figure [5] shows ^([/, n, 6) 
for 9 = 1,2, where U is defined as in Proposition 14.7 1 with e = 0.5. Note that the numerical computation 
actually shows that 

^{U,n,9) = 0{n), (5.9) 

which is indeed a very pleasant surprise. In fact, due to the "stair case" growth shown in Figure |5] the 
growth is actually better than what ( 15.91 ) suggests. The question is whether this is a particular quality of 
the Haar wavelet, or that one can expect similar behavior of other types of wavelets. The answer to this 
question will be the topic of future papers. 

Note that Figure |5] is interpreted as follows: provided m > 4.9n, for example, we can expect this 
method to reconstruct g to within an eiTor of size (1 + 0)||P^/3||, where 6* = 1 in this case. In other 
words, the error is only two times greater than the best approximation to g from the finite-dimensional 
space consisting of the first n Haar wavelets. 

Having described how to determine conditions which guarantee existence of a reconstruction, in the 
next section we apply this approach to a number of example problems. First, however, it is instructive to 
confirm that these conditions do indeed guarantee stability of the recontruction procedure. In Figure [6] we 
plot ||(e^)~^|| against n (for e = 0.5), where A is formed via (14.191 ) using Haar wavelets with parameter 
m = [4.9n]. As we observe, the quantity remains bounded, indicating stability. Note the stark contrast 
to the severe instability documented in Figure |2] 



6 Examples 

In this final section, we consider the application of the generalized sampUng theorem to several examples. 

6.1 Reconstruction from the Fourier Transform 

In this example we consider the following problem. Let / G L^ (M) be such that 

/ = Tg, supp(fi() C [-T,T]. 
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50 100 150 200 250 300 350 
Figure 6: The quantity ||(e^)~^|| against n = 2, 4, . . . , 360. 



We assume that we can access point samples of /, however, it is not / that is of interest to us, but rather g. 
This is a common problem in applications, in particular- MRI. The NS Sampling Theorem assures us that 
we can recover g from point samples of / as follows: 



oo 

E 

i=— oo 



fine) e 



2mne- 



1 

2r' 



where the series converges in L? norm. Note that the speed of convergence depends on how well g can be 
approximated by the functions e^'^'"'^', n G Z. Suppose now that we consider the function 

g[t) = cos(27rt)x[o.5,i](*)- 
In this case, due to the discontinuity, forming 

AT 



9N 



E /(^ 



e 2^ J\'ne)e 

n=~N 



27rme- 



e = \, iVGN, 



(6.1) 



may be less than ideal, since the convergence gj\f ^ g as N ^ oo may be slow. 

This is, of course, not an issue if we can access all the samples {f{ne)}nez, however, such an as- 
sumption may be overly optimistic in application. Moreover, even if we had access to all samples, we are 
limited by both processing power and storage to taking only a finite number. 

Suppose that we have a more realistic scenario: namely, we are given the finite collection of samples 



7?/ = {f{-Ne)J{{-N + l)e), . . . , f{{N - l)e)J{Ne)}, 



(6.2) 



with N = 900 and e = |- The task is now as follows: construct the best possible approximation to g 
based on the vector rjf. We can naturally form g^ as in (16.11 ). This approximation can be visualized in the 
diagrams in Figure |2l Note the rather unpleasant Gibbs oscillations that occur, as discussed previously. 
The problem is simply that the set {e^'^'"'^}„gz is not a good basis to express g in. Another basis to use 
may be the Haar wavelets {ipj} on [0, 1] (we do not claim that this is the optimal basis, but at least one 
that may better capture the discontinuity of g). In particular, we may express g as 



oo 



z^/^jV'j, 



/3 = {/3i,/32,...}G/'(N). 
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Figure 7: The upper figures siiow (77V (left), gn,m (middle) and g (right) on the interval [0, 1]. The lower 
figures show g^ (left), gn,m (middle) and g (right) on the interval [0.47, 0.57]. 



We will now use the technique suggested in Theorem 14.81 to construct a better approximation to g based 
on exactly the same input information: namely, rjj in (16.21 ). Let U be defined as in (14.181 ) with e = 1/2 
and let n = 500 and m = 1801. In this case 



PnU*PmUPn\p„H 



PnU*PmUPn\pM] PnU* P, 



< 0.6169, 

< 0.7854 



Define /? = {/5i, . . . , /?„} by equation (14.191) . and let gn,m = Yll=i Pji^j- The function gn,m is visualized 
in Figure |7] Although, the construction of gj\i and gn,m required exactly the same amount of samples of 
/, it is clear from Figure |7] that gn,m is favorable. In particular, approximating g by gn,m gives roughly 
four digits of accuracy. Moreover, had both n and m been increased, this value would have decreased. In 
contrast, the approximation g^ does not converge uniformly to r^ on [0, 1]. 

6.2 Reconstruction from Point Samples 

In this example we consider the following problem. Let / G L^ (M) such that 

K 

f = ^g, g{x) = ^aj>i(a;) +sin(27rx)x[o.3,o.6](a;), 
j=i 

for K = 400, where {V'il are Haar wavelets on [0,1], and {aj}^^^ are some arbitrarily chosen real 
coefficients in [0, 10]. A section of the graph of / is displayed in Figured! The NS Sampling Theorem 
yields that 



fit) = E / (^) sinc(2t 



k), 



fc=— 00 



24 




20 

10 



-10 

-20 




■■''l'' - l'l'l' 4 l n ttH ton || ntn|r 



-5000 5000 -5000 5000 

Figure 8: The figure siiows Re(/) (left) and Im(/) (right) on the interval [—5000, 5000]. 

where the series converges uniformly. Suppose that we can access the following pointwise samples of /: 

r^f = {f{-Ne)J{{-N + l)e), . . . , f{{N - l)e)J{Ne)}, 

with e = ^ and A^ = 600. The task is to reconstruct an approximation to / from the samples rjf in the 
best possible way. We may of course form 

N 



" / 7 \ 

fN{t)= Yl fl^)smc{2t-k), N = 600. 



k=-N 

However, as Figure |9] shows, this approximation is clearly less than ideal as f{t) is approximated poorly 
for large t. It is therefore tempting to try the reconstruction based on Theorem 14. 8 1 and the Haar wavelets 
on [0, 1] (one may of course ti"y a different basis). In particular-, let 

n 

/ = 5Z^i-^^i' ^ = 500, 



where 



A/3 = PnU*PmVf, A = PnU*PmUPn\p„n, 



with m = 2N + 1 = 1201 and U is defined in (14.181 ) with e = 1/2. A section of the errors \f — fj^\ and 
1/ — /I is shown in Figure |9] In this case we have 



PnU*PmUPn\pM 



PnU*Pm.UPn\p„H] PnU* P, 



-1 



< 0.9022, 

< 0.9498. 



In particular, the reconstruction / is very stable. Figure |9] displays how our alternative reconstruction is 
favorable especially for large t. Note that with the same amount of sampling information the improvement 
is roughly by a factor of ten thousand. 

7 Concluding Remarks 

The framework presented in this paper has been studied via the examples of Haar wavelets and Legendre 
polynomials. Whilst the general theory is now well developed, there remain many questions to answer 
within these examples. In particular. 
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Figure 9: The figure shows the en^or \f — JnI (left) and \f — f\ (right) on the interval [—5000, 5000]. 



(i) What is the required scaling of m (in comparison to n) when the reconstruction basis consists of 
Legendre polynomials, and how well does the resulting method compare with more well-established 
approaches for overcoming the Gibbs phenomenon in Fourier series? Whilst there have been some 
previous investigations into this particular approach |[l6l[T9l, we feel that the framework presented 
in this paper, in particular the estimates proved in Theorem 14.11 are well suited for understanding 
this problem. We are cun^ently investigating this possibility, and will present our results in a future 
paper. 

(ii) Whilst Haar wavelets have formed been the principal example in this paper, there is no need to 
restrict to this case. Indeed, Theorem 15.11 provides a first insight into using more sophisticated 
wavelet bases for reconstruction. Haar wavelets ai^e extremely simple to work with, however the 
use of other wavelets presents a number of issues. In particular, it is first necessary to devise a 
means to compute the entries of the matrix [/ in a more general setting. 

In addition, within the case of the Haar wavelet, there remains at least one open problem. The 
computations in Section [5TT] suggest that n i— ;• ^{U,n,6) is bounded by a linear function in this 
case, meaning that Theorem 15. II is overly pessimistic. This must be proven. Moreover, it remains 
to be seen whether a similar phenomenon holds for other wavelet bases. 

(iii) The theory in this paper has concentrated on linear reconstruction techniques with full sampling. A 
natural question is whether one can apply non-linear techniques from compressed sensing to allow 
for subsampling. Note that, due to the infinite dimensionality of the problems considered here, the 
standard finite-dimensional techniques are not sufficient. 
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